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ABSTRACT 

Using fully three-dimensional hydrodynamic simulations, we investigate the effect of the Coriolis 
force on the hydrodynamic and observable properties of colliding wind binary systems. To make the 
calculations tractable, we assume adiabatic, constant velocity winds. The neglect of radiative driving, 
gravitational deceleration, and cooling limit the application of our models to real systems. However, 
these assumptions allow us to isolate the effect of the Coriolis force, and by simplifying the calculations, 
allow us to use a higher resolution (up to 640^) and to conduct a larger survey of parameter space. 
We study the dynamics of coUidng winds with equal mass loss rates and velocities emanating from 
equal mass stars on circular orbits, with a range of values for the ratio of the wind to orbital velocity. 
We also study the dynamics of winds from stars on elliptical orbits and with unequal strength winds. 
Orbital motion of the stars sweeps the shocked wind gas into an Archimedean spiral, with asymmetric 
shock strengths and therefore unequal postshock temperatures and densities in the leading and trailing 
edges of the spiral. We observe the Kelvin-Helmholtz instability at the contact surface between the 
shocked winds in systems with orbital motion even when the winds are identical. The change in 
shock strengths caused by orbital motion increases the volume of X-ray emitting post-shock gas with 
T > 0.59 keV by 63% for a typical system as the ratio of wind velocity to orbital velocity decreases to 
Vw/Vo — 2.5. This causes increased free-free emission from systems with shorter orbital periods and 
an altered time-dependence of the wind attenuation. We comment on the importance of the effects of 
orbital motion on the observable properties of colliding wind binaries. 

Subject headings: binaries: general — hydrodynamics — stars: early-type — stars: winds, outflows 
— stars: Wolf-Rayet — X-rays: stars 



1. INTRODUCTION 

A large proportion (at least 39%) of Wolf-Rayet (WR) 
stars in the solar neighborhood have been observationally 
confirmed to lie in binary systems (van der Hucht 2001). 
These stars have dense, highly-supersonic, radiatively- 
driven winds that, upon colliding with the wind from 
a companion, produce strong shocks that compress and 
heat the gas to temperatures high enough to produce 
substantial X-ray flux. These colliding wind binaries 
(CWBs), which can contain a pair of WR stars, a WR 
star with an OB supergiant companion, or a pair of OB 
stars, were first proposed by Prilutskii & Usov (1976) 
and Cherepashchuk (1976). This X-ray excess has since 
been observationally confirmed (Pollock 1987). 

For more than a decade, the hydrodynamics of CWBs 
have been investigated using two-dimensional (axisym- 
mertric) hydrodynamical simulations (e.g. Luo, McCray, 
& Mac Low 1990; Stevens, Blondin, & Pollock 1992) in 
order to compute synthetic X-ray spectra for comparison 
to observations. A variety of physical effects that are im- 
portant for the dynamics of the winds have been iden- 
tified, including radiative inhibition (Stevens & Pollock 
1994) and sudden radiative braking (Gayley, Owocki, & 
Cranmer 1997). The former is the effect of the radia- 
tion field of one star decreasing the radiative accelera- 
tion of the wind from the other, while the latter is where 
the radiation field of one star is intense enough to ac- 
tually cause a net deceleration of the other star's wind. 
While including such effects leads to more realistic mod- 
els of CWBs (e.g. Henley, Stevens, & Pittard 2005), most 
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current simulations are typically still done in only two- 
dimensions, assuming axisymmetry around the line join- 
ing the two stars. This requires neglecting the Coriolis 
force resulting from orbital motion. Including the Corio- 
lis force requires fully three-dimensional hydrodynamical 
models which, up to now, have been largely untractable. 

In this paper, we investigate the effect of orbital mo- 
tion on the hydrodynamic and observable properties 
of the wind collision region by performing fully three- 
dimensional hydrodynamic simulations of CWBs. We 
conduct a parameter survey in orbital velocity for stars 
on circular orbits, and also present models for stars on 
elliptical orbits and with unequal strength winds. 

Our focus in this paper is the three-dimensional dy- 
namics of the wind collision region, thus to simplify the 
problem we neglect a number of physical processes that 
can be important in real systems, but would serve to 
complicate the interpretation of our results. In particu- 
lar, we neglect radiative cooling of the postshock gas and 
do not try to model the radiative driving of the winds. 
Instead of explicitly including gravity in our simulations, 
we assume that gravitational deceleration balances ra- 
diative acceleration in the unshocked winds and give the 
winds a constant velocity profile. This velocity is chosen 
to be the unshocked wind velocity given by a beta veloc- 
ity law at the distance of the stagnation point (see eq. 
[5] for more details). We discuss in [2] possible implica- 
tions of neglecting gravitational force on the unshocked 
wind of the companion star as well as the post-shock 
gas. Since we do not model radiative effects, there will 
be no sudden radiative braking or radiative inhibition in 
our simulations. Sudden radiative braking is thought to 
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occur only with strongly unequal winds (Gayley et al. 
1997), which we will not consider here. 

While for most early-type binaries it is important to 
consider radiative cooling of the postshock gas by line 
emission, we leave this to future research. We have cho- 
sen a sample system that could be composed of a pair 
of 0-stars, which we find is neither purely adiabatic nor 
strongly radiative (see using the appropriate for- 

mula in Antokhin, Owocki, & Brown (2004). We will 
comment on how including more realistic radiative ef- 
fects on the winds might change our results. 

It is important to note ours are not the first fully three- 
dimensional simulations of CWBs; Walder (1998) showed 
that orbital motion can produce Archimedian spirals sim- 
ilar to those visible in infrared images of WR 98a (Mon- 
nier, Tuthill, & Danchi 1999) and WR 104 (Tuthill, Mon- 
nier, & Danchi 1999). 

The organization of this paper is as follows. In ij2] we 
describe the details of our simulations and the diagnostics 
we use to interpret our results. In 321 we present anal- 
ysis of the post-shock temperatures, emissivity, column 
density, and instabilities that result for colliding winds 
from identical stars on circular orbits, with a range of 
values for the ratio of wind velocity to orbital velocity. 
We also consider systems with eccentric orbits and un- 
equal winds. Finally, in iJH we discuss the implications 
by applying our results to WR 20a. 

2. METHOD 

We present three-dimensional hydrodynamic simula- 
tions conducted on a cartesian grid with the Athena code 
(Gardiner & Stone 2005, 2006). We solve the equations 
of hydrodynamics 



dp 

dt 



V • (pv) = 0, 



^ -f V • (pvv + P) = 0, 



and 



dE 
'dt 



V-[iE + P)v] = 0, 



where E is the total energy density, 
E = e + —pw, 



(1) 
(2) 

(3) 
(4) 



e is the internal energy density, and we adopt an ideal 
gas equation of state, P = (7 — l)e, where 7 = 5/3 is the 
ratio of specific heats. We use the Roe Riemann solver 
augmented with the H-correction of Sanders, Morano & 
Druguet (1998) to prevent the carbuncle phenomenon 
(Quirk 1994). We use outflow boundary conditions and 
run our simulations for up to three orbital periods. Most 
of our simulations were run on 256^ grids with boxes of 
length L — 2.5a AU (where the seperation between the 
stars is a AU) , however we also ran three simulations on 
640^ grids, two of them with L — 6.25a AU and the third 
with L = 5a AU. We will define the parameter a in the 
next section. 

2.1. Initial Conditions 

Our stars are given a steady, spherically-symmetric, 
radial wind by forcing the correct density, pressure, and 



TABLE 1 

Simulation Parameters 



Parameter 




Scaling Used 


Semi-major axis of system 


a = 


a/1 AU 


Wind velocity 


V = 


t)»/1000 kms-i 


Mass loss rate 


V = 


M/lO-^Mgyr-l 


Wind mach number 


Mo 


= 30 



velocity profiles onto a small region of the grid. The 
density and pressure of the wind are given by 



p(r) =po(^) 



and 



PW = P.(^) 



10/3 



(5) 



(6) 



Our wind is generated with a constant velocity, Vw , taken 
to be the velocity at the stagnation point in the circular- 
orbit case. The constants po and Po are defined to give 
the desired mass loss rate, M, and wind mach number, 
A4q, at the edge of the mask, tq, using 



M = 4:nrlpoVn 



and 



Mo 



Po 

iPo'' 



(7) 



(8) 



where ro is some fiducial radius, taken to be ro — 
0.195aAU. The density, pressure, and velocity in the 
simulations are in an arbitrary system of units. 

The wind and orbital parameters in our simulations 
can been parameterized by the variables a, v, and 77, de- 
fined in Table [TJ These quanities, which are generally of 
order unity for early- type binary systems, can be varied 
in order to apply our results to a variety of systems. 

Initially, the entire computational domain is filled 
with a stationary ambient medium with density Pamb = 
1.09 X 10~^''77a~^:^~-^ gcm"'^ and sound speed Cs.amb = 



21.0j/kms 



The masks are then initialized with a ra- 



dius 2ro. The properties of the ambient medium were 
chosen for numerical, not physical, reasons. We evolve 
the system for much longer than a wind-crossing time 
to make sure the ambient medium has been completely 
driven off the grid before starting our analysis. See the 
appendix for the details of how we generate our wind. 

2.2. Systems Investigated 

To investigate the effect of orbital motion on both the 
hydrodynamic and observable properties of the collision 
region, we take a binary system with a circular orbit and 
vary the ratio of the wind velocity to orbital velocity. 
The mass of the stars and the orbital period are varied 
so as to keep the semi-major axis of the orbit constant 
while varying the orbital velocity. 

The sample system to which we will apply our results 
has a = 0.69, ly = 0.85, and = 0.065. This gives a = 
0.69 AU, M = 6.5 x lO-^Mgyr"!, and = 850km/s 
at the stagnation point. Given a beta velocity law. 



v{r) 



(9) 
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TABLE 2 
Circular-orbit 
Simulation 
Parameters 



Name 




5 


oo 


CIO 


10.0 


C5 


5.0 


C3.5 


3.5 


C2.5 


2.5 



with /3 = 1, we find Vx, ~ 1160 kms~ for its terminal 
velocity, assuming i?, — 20Rq. This set of parameters 
gives for each of the stellar masses 

= 2.25 X 10^ aiy^(^) ^ Mg (10) 

= 1.12x10^(^)"'mo. (11) 

It should be noted that stars with masses above those 
given by Vw/Vo « 3.6 for this sample system have not 
yet been detected. Using the method of Antokhin et al 
(2004) to determine if a system with these parameters is 
radiative (a < Urad) or adiabatic {a > arad), we find 

arad = 4.66^ = 0.683, (12) 

meaning that our system is very slightly on the adiabatic 
side. Despite neglecting radiative effects in our simula- 
tions, we feel there is still something to be learned from 
studying this hypothetical set of systems. 

We start by presenting the results of circular-orbit sim- 
ulations with Vw/Vo = oo and Vw/Vo = 2.5 in a large 
box {L ~ 6.25q!AU) on a 640'^ grid. We then compare 
a set of simulations with varying Vw/Vo, whose names 
and velocity ratios are given in Table [51 in a smaller box 
(L = 2.5q!AU) on a 256^ grid- This preserves cell size 
between the large- and small-box runs. 

We compare two additional simulations with eccentric 
orbits to our corresponding circular-orbit simulations. 
Simulations SE2.5 and E2.5, respectively, have identi- 
cal properties to simulations S and C2.5, except with 
eccentricity e — 0.2. In simulation SE2.5, the separation 
of the stars is varied over time with zero orbital veloc- 
ity perpendicular to the line of centers, as was done by 
Pittard (1998). 

Finally, we consider simulations similar to C2.5 but 
with unequal winds. Simulation CM2.5 has identical 
parameters to C2.5, except that, for the secondary (on 
the negative x-axis when the simulation is initialized), 
?72 — (2/3)77, giving a lower mass-loss rate. Simulation 
CW2.5 also has identical parameters to C2.5, except 
that, for the secondary, 1/2 = (2/3)i^, giving a lower wind 
velocity. The modified wind velocity and mass-loss rate 
were chosen so as to give a wind momentum ratio of 1.5 
in all unequal wind simulations. 

2.3. Diagnostics 

We analyze our simulations by comparing post-shock 
temperatures, free-free emission from post-shock gas 
above a minimum temperature, and column density to 
key locations in the system. So that our results can be 
scaled to other systems, we leave most of our results pa- 
rameterized in terms of a, v, and 77, described in 



as well as fl, the mean particle mass in amu, and Z, the 
charge of the ion, but substitute for Mo- Stevens et al. 
(1992) give fl = 0.6 for solar abundances, fi — 1.3 for WN 
stars, and /2 = 1.4 for WC stars, assuming full ionization. 

Due to the scalings chosen for our simulations, one can 
convert the densities from our simulations to physical 
units using 

p = 1.09 X lO-i^a-^T^-Vgcm-^, (13) 
and the temperatures by 

T = 0.0106/21/2 CkeV, (14) 
P 

where P and p are the pressure and density, respectively, 
in simulation units. 

We calculate the power radiated per unit volume by 
free- free emission over all frequencies and directions using 

Ay/ - 2.1 X 10-«Z2/2-2r^/^p2_^^ergs-icm-3, (15) 

where p-ie is the density in units of lO^^^gcm"'^. We 
have assumed the frequency-averaged Gaunt factor to 
be (gff) = 1.2 with the understanding that this will 
only give an accuracy of about 20% (Rybicki & Light- 
man, 1979). We integrate this equation over the post- 
shock volume to find the total power radiated by free- 
free emission, Vff. The kinetic luminosity of each wind 
is Lw = 3.15 X lO^'^rjv^ ergs'i. 

To see how the light curve would be attenuated by the 
wind, we calculate the column density to each star and to 
the center of mass of the system, the source of the hardest 
X-ray emission from the post-shock gas. We use lines of 
sight inclined 25° from the normal to the orbital plane, 
whose projections onto the orbital plane are aligned with 
the zero-phase line of centers. Integrating from the edge 
of the box to either the center of mass or the edge of the 
stellar mask, we find for the column density 

N = 8.95 X lO^V"^ J P-wdvAU cm"^- (16) 

Since we are only integrating from the edge of the box, we 
will miss the contribution to the column density from the 
edge of the box to the observer. The contribution from a 
spherically-symmetric, unshocked wind outside the box, 
along the line of sight to one of the stars, would be 

AN = 3.64 X lO^^^LXu cm-^ (17) 

For our sample system in a small box (L = 2.5a AU), 
this gives AN = 1.62 x lO^^/I"^ cm"^. In the largest box 
{L = 6.25a AU), we find AiV = 6.46 x lO^^/i^i cm^^. 

3. RESULTS 

3.1. Identical Winds on Circular Orbits 

Figure [T] shows the density for the large-box {L = 
6.25a AU) runs of S and C2.5 in three orthogonal cut 
planes through the center of mass of the system. Figure 
[2] shows the temperature in the same three planes. The 
shape of the shock front for C2.5, with Vw/Vo — 2.5, dif- 
fers substantially from that for the stationary star case, 
S. The shocks wrap around the stars, breaking the ax- 
isymmctry that is present when the stars are held sta- 
tionary. The absence of shocked material in the two cir- 
cular regions surrounding the stars in the slices normal 
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Fig. 1. — Density in large-box simulation S (a) in orbital plane, 
(b) through the stars normal to the orbital plane, and (c) normal 
to the previous two slices. Density in large-box simulation C2.5 
(d) in orbital plane, (e) through the stars normal to the orbital 
plane, and (/) normal to the previous two slices. The color scale 
is logarithmic, with a^urj-ip-iQ from 1 (white) to 700 (red). The 
tick marks are spaced aAU. 
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Fig. 2. — Temperature in large-box simulation S (a) in the 
orbital plane, (b) through the stars normal to the orbital plane, 
and (c) normal to the previous two slices. Temperature in large- 
box simulation C2.5 (d) in the orbital plane, (e) through the stars 
normal to the orbital plane, and (/) normal to the previous two 
slices. The color scale is logarithmic, with p,~^u~^T]^c\/ from 0.05 
(white) to 2.5 (red). The tick marks are spaced aAU. 

to the orbital plane is a result of the curvature of the 
shock. The shocked material that would be projected 
onto those regions is on the near side of the left star and 
on the far side of the right star, relative to the viewer, 
as is apparent from the slice in the orbital plane. The 
three-dimensional structure of the post-shock region is 
more visible in Figure [3l which traces the shock fronts 
using a surface of constant temperature. Small spheres 
mark the size and position of the stellar masks. There is 
63% more volume at T > OMfliy'^ keV in C2.5 than in S 
but only 20% more at T > 0.875/2;^^ keV. 

The maximum post-shock density and temperature are 
not significantly different between C2.5 and S when the 
magnitude of the time-variation in C2.5 is considered. 
The temperatures on the leading side of the shock in 
C2.5 are higher than on the trailing side; conversely, the 
densities are lower. The thickness of the post-shock re- 
gion between the stars is the same for both C2.5 and S. 
In both cases, the post-shock region is centered on the 
center of mass of the system. The numerical methods 
used in the Athena code keep the shocks very thin, lead- 
ing to a stair-stepping of the density and temperature as 
the curved shock crosses the Cartesian grid. This stair- 
stepping is clearly visible in Figure IH a slice through the 
orbital plane of the temperature in the small-box simu- 
lation C5. 




Fig. 3. — Shock structure for large-box circular-orbit simulation 
with Vw/Vo = 2.5. The surface shown is T = 0mi25p.u^ keV. 
Small spheres mark the size and location of the stellar masks. The 
stars are moving counter-clockwise in the x-y plane. 

While there is no evidence of the Kelvin-Helmholtz (K- 
H) instability at the contact surface bewteen the shocked 
winds in S (where the stars are held stationary and the 
wind velocities are identical), K-H instability is clearly 
visible, as in Figure 31 in simulations where the stars 
are orbiting, even though the winds are identical. Al- 
though over most of the stellar orbits the effect of the 
K-H instability on the postshock gas is very small, in- 
termittently large rolls form which can affect observable 
quantities. For example, the power emitted from gas 
with T > 0.875^1^^ keV varies by ~ 2% as a large roll 
is advected off the grid in the small-box simulation C5. 
The instability is seeded by grid noise associated with 
the representation of thin shocks on the Cartesian grid 
(the stair-stepping observed in Figure [3]). The amplitude 
of the grid noise is correlated with the phase of the orbit 
(direction of the shock with respect to the grid), there- 
fore we see periodicity in the K-H rolls at one-quarter 
of an orbit. Although the K-H rolls in our calculations 
are seeded by grid noise, the effect is real, and could be 
seeded by variability or dumpiness in the winds. The 
stair-stepping also introduces grid noise into the volume 
occuppied by gas at various temperatures. For this rea- 
son, we try to compare the results for rapidly orbiting 
stars to model CIO rather than model S as much as pos- 
sible in our parameter survey, since in CIO the motion of 
the shocks over the grid averages out the grid noise, yet 
the orbital motion is small enough to not significantly 
affect the structure compared to the stationary case S. 

The power emitted by free-free emission from 
gas with T > 0.63/2j/2keV is Vff = 1.9 x 
1036z2/2~3/2^2Q,~ij,-igj.gg-i (J2.5. This gives the 
ratio of the free-free power to the wind luminosity 
0.29Z^/i~'^/^77^a~^^~^. For our binary 0-star system, 
this evaluates to 2.1 x 10^'^Z^/2~'^/^, which supports not 
incorporating cooling into our simulations for this par- 
ticular system. The time-average of the power radiated 
by free-free emission from gas with T > 0.63/ii^^ keV is 
13% higher for C2.5 than for 5*. The standard deviation 
is 0.24% of the mean for the power emitted from C2.5. 
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Fig. 4. — Temperature in the orbital plane of small-box sim- 
ulation C5. The image has been enhanced to make the effect of 
the Kelvin-Helmholtz instability on the contact surface more easily 
visible. The tick marks are spaced 0.5a AU. 

Figure [5] shows volume as a function of post-shock tem- 
perature averaged over one orbital period. Post-shock 
temperature decreases rapidly with distance from the line 
of centers, leaving most of the post-shock volume at low 
temperatures. Error bars, resulting from both grid ef- 
fects and the K-H instability, are shown but are barely 
visible except at the lowest temperatures plotted. As in- 
stabilities are seeded near the center of mass of the sys- 
tem, the volume at the highest temperatures is affected 
first. A wave with an amplitude of a few percent passes 
through the temperature distribution, moving from high 
to low temperatures, as the gas moves away from the line 
of centers and off the grid. 

As expected, the highest post-shock temperatures are 
near the line of centers. Using the density and pressure 
profiles described in §2.11 we derive an unshocked tem- 
perature and mach number at the stagnation point in 
the stationary-star simulation that would yield a maxi- 
mum post-shock temperature of T = keV. This 
temperature is marked by the long dashed vertical line 
in Figure [HI The high-temperature tail of the distri- 
bution extends beyond this predicted maximum, likely 
due to multi-dimensional effects not taken into account 
in our estimate. To prevent our results from having 
a dependence on orbital phase due solely to the cubi- 
cal shape of the computational domain, we choose to 
consider in our analysis only cells with post-shock tem- 
peratures T > keV for the large-box runs and 
T > 0. 875 fLi/^ keV (marked with the short dashed verti- 
cal line in Figure [5]) for the small-box runs. 

Figure [6] shows volume as a function of power 
per unit volume due to free-free emission from gas 
with T > keV. In C2.5, gas with A// > 

4.5 X lO^^Z^/i^^/^Ty^a^**!/^^ ergs"^ cm~^, all has T > 
0.875/2:/^ keV, therefore using that cutoff instead of T > 
keV will not cause us to miss the high-emissivity 
gas. While most of the post-shock volume at low temper- 
atures emits little power per unit volume, the integrated 
power from this volume, as shown by Figure [71 can still 
contribute significantly to the total emission. The power 



Fig. 5. — Histogram of volume within post-shock temperature 
bins of width AT = Q.0212flu^ keV, averaged over one orbit, from 
large-box runs of C2.5 (solid) and S (dotted). The vertical short 
dash line marks T = keV, the temperature cutoff that 

will be used for the circular-orbit parameter survey analysis. The 
vertical long dash line marks T = LdGfiu^ keV, the expected max- 
imum temperature calculated analytically. 
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Histogram of volume within bins of width AA 
-3/2^2^- 



1.25 X 10~ Z r; a~ ergs" cm~^ in volume-averaged 

power per unit volume due to free-free emission from gas with 
T > keV for large-box C2.5 (solid), averaged over one or- 

bit, and large-box S (dotted). The peak power per volume is higher 
when the stars are orbiting than when they are stationary. 

emitted from gas with T > keV is only 79% 

of that emitted from all gas with T > O.Q3p,v^ keV in 
the large-box C2.5. Similarly, the power emitted from 
gas with T > 0.875^1/^ keV is only 82% of that emitted 
from all gas with T > 0.63fliy^ keV in the large-box S. 
The quantitative change between observables measured 
will depend on the cutoff used, but since we are using 
the same cutoff for each simulation when making direct 
comparisons, we will still get a qualitative feel for the 
effect as well as an esimate of its magnitude. 

Column density, as described in ij2.31 is plotted in Fig- 
ure [H The column density to the center of mass is deter- 
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Fig. 7. — Histogram of power emitted from bins of width 
AAff = 1.25 X 10~^Z'^p.~^/'^ri^a~^L'~^ evgs~^ cm~^ in volume- 
averaged power per unit volume due to free-free emission from gas 
with T > 0.63p,u^ keV for large-box C2.5 (solid), averaged over one 
orbit, and large-box S (dotted). There is more integrated emission 
from low power per volume but less from high power per volume 
when the stars are stationary compared to when they are orbiting. 

mined mostly by the density of the post-shock material 
and the thickness of the post-shock region along the line 
of sight. At phases 0.5 and 1.0, the line of sight falls in 
the cut plane shown in Figure [TJb, at an angle 25° from 
vertical. At phases 0.25 and 0.75, the line of sight falls 
in the cut plane shown in Figure [If. The column density 
at phases 0.25 and 0.75 is higher than at 0.5 and 1.0 be- 
cause, although the line of sight passes through slightly 
less post-shock material, it passes through more post- 
shock material close to the center of mass of the system, 
where the mass density is highest. When a smaller box 
is used, the column density to the center of mass changes 
most near the minima. The post-shock gas that was near 
the edge of the box in the large-box runs is not visible in 
the small-box runs. Even though the values change, the 
qualitative features remain. 

The column density to the primary and secondary 
when the orbit is circular and the winds are identical 
is essentially the same but with a phase shift of 0.5. It 
is highest near when the stars are on the far side of the 
center of mass relative to the observer, near phase 0.0 
for the primary and 0.5 for the secondary. When the 
stars are equally distant from the observer, the column 
density is higher to the star that is moving toward the 
observer, the primary at phase 0.25 and the secondary 
at phase 0.75, because the post-shock material leads the 
star. When a smaller box is used, the column density 
to the stars is affected at all phases. Instead of contin- 
uously decreasing after one peak until the next peak is 
approached, in a small box the column density flattens 
out well before the next peak. 

3.2. Circular-orbit Parameter Survey 

Figure [5] shows density contours both in the orbital 
plane and normal to it, as well as temperature contours 
in the orbital plane, for five different simulations that 
survey the ratio of wind to orbital velocity. All use stars 
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Fig. 8. — Column density on the line of sight described in i|2.3l to 
the center of mass (top), primary star (middle), and secondary star 
(bottom), for C'2.5 in a large box (solid) and small box (dotted). 
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Fig. 9. — Logarithmic density contours in the orbital plane (top) 
and normal to the orbital plane (middle) and linear temperature 
contours in the orbital plane (bottom) for (left to right) small-box 
circular-orbit simulations 5, CIO, C5, C3.5, and C2.5. The smaller 
size of the computational domain (L = 2.5oAU) means the spiral 
pattern seen at larger distances from the stars is not captured. 

on circular orbits in a box of size L = 2.5a AU. The 
stars are moving counter-clockwise, as before. The shock 
fronts curve around the stars more and more as the or- 
bital period decreases, both in the orbital plane and nor- 
mal to it. The smaller box crops off the spiral patterns 
visible in the large-box simulations (Figures [1] and [2]) ■ 
The maximum post-shock density shows a nearly mono- 
tonic decrease of ~ 6% as the velocity ratio decreases to 
Vw/Vo = 2.5, while the maximum post-shock tempera- 
ture doesn't change significantly. 

The power radiated by free- free emission from gas with 
T > 0.875/1^^ keV increases as orbital period decreases. 
The time-average of the power emitted from CIO is only 
0.085% higher than from S. This justifies using CIO as 
a reference model instead of S in the discussion to follow 
since S is stationary on the grid and will have more grid 
noise. The time-average of the power emitted from C5 
is 1.9% higher, from C3.5 is 4.0% higher, and from C2.5 
is 7.7% higher than from CIO. 

Figure [To] shows the fractional difference in post-shock 
volume, averaged over one orbital period, as a function 
of post-shock temperature for the circular-orbit simula- 
tions when compared to simulation ClO, with the low- 
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Fig. 10. — Histogram of fractional difference in volume, averaged 
over one orbit, within post-shock temperature bins of width AT = 
0.0212/i!^^ keV for simulations C5 (short dash), C3.5 (long dash), 
and C2.5 (dash dot) relative to CIO (dotted). There is more volume 
at lower temperatures and less at higher temperatures when the 
orbital period is shorter. 



est non-zero orbital velocity. For temperatures below 
T K. keV, the volume is larger for simulations 

with shorter orbital periods because the post-shock tem- 
peratures are higher on the leading side of the post-shock 
region than they would be if the stars were stationary. 
For higher temperatures, the volume is smaller for sim- 
ulations with shorter orbital periods because the wind 
impacts the shock front along the line of centers at a 
more oblique angle. 

We histogram the volume in bins of width AA/y = 
1.25 X \Q^^Z^\i~^l'^"(f'a~^v'^ ergs^^ cm^'^ in power per 
unit volume due to free-free emission from gas with T > 
0.875/i:/2 for C5, C3.5, and C2.5 compared to CIO. 
CIO has the highest maximum power per unit volume, 
although C3.5 and C2.5 have more volume near that 
value. At intermediate values of power/volume, in the 
range 1.5 x 10"^ < f^^g.Z-'^lt'l'^jf'^a^v < 4.5 x lO'^, 
the volume in a given bin is an average of 2.7% higher 
for C5, 6.8% higher for C3.5, and 15% higher for C2.5 
than for CIO. The power emitted from the trailing side 
of the post-shock region is higher than from the leading 
side when there is orbital motion. 

The ratio of wind velocity to orbital velocity impacts 
the column density to the hot post-shock gas as well. 
Figure [TT] shows column density for each of the circular- 
orbit simulations. The column density to the center of 
mass of the system, where there will be the hardest X-ray 
emission from the post-shock gas, peaks twice per orbital 
period. While the column density should peak when the 
stars are half way between conjunctions when there is no 
curvature of the shocked gas due to the Coriolis force, 
as the orbital period decreases and the Coriolis force be- 
comes stronger, the peak column densities are delayed 
and the peaks become less symmetrical, increasing more 
slowly and then decreasing more rapidly. The peak is de- 
layed by a phase of 0.07 for C2.5 compared to CIO and 
the difference between the peak and minimum column 
density increases by 13% between CIO and C2.5. The 
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Fig. 11. — Column density on the line of sight described in i|2.3l 
to the center of mass (top), primary star (middle), and secondary 
star (bottom). 

FWHM remains roughly constant as the orbital velocity 
is varied. 

When the orbital period is very long, the column den- 
sity to the stars is nearly symmetrical in time, but as the 
orbital period decreases and the curvature of the shocked 
region becomes stronger, the column density becomes 
very asymmetrical. As the orbital period decreases, the 
column density increases more rapidly leading up to the 
peak and then decreases less rapidly after the peak, caus- 
ing an increase in the FWHM of the peaks of 31% from 
CIO to C2.5. The column density of the star that is fa- 
ther away from the observer increases quickly as the stars 
approach conjunction. The peak column density should 
occur at conjunction when the shock fronts aren't curved 
by the Coriolis force, but when they are curved, the peak 
occurs after conjunction. For CIO, the peak occurs at a 
phase only 0.02 past conjunction, but for C2.5, the peak 
is delayed to a phase 0.07 past conjunction. The dif- 
ference between the peak and minimum column density 
increases by 22% between CIO and C2.5. The fraction 
of the orbit during which the column density is signif- 
icantly above its minimum value is higher for systems 
with shorter orbital periods. It is only 41% for CIO but 
increases to 57% for C2.5. 

3.3. Unequal Winds on Elliptical Orbits 

To show how a combination of an elliptical orbit and 
unequal winds affects the shock structure, we present 
images of a large-box [L — 5a AU) elliptical-orbit simu- 
lation with unequal mass-loss rates. The mass-loss rate 
of the primary has been increased by a factor of 1.5 com- 
pared to the secondary and the eccentricity of the orbit 
is e = 0.2. 

Figure [12] shows the density in three orthogonal cut 
planes that pass through the center of mass of the system. 
These planes are shown at three different phases of the 
orbit. Figure [T3] shows the temperature in the same three 
planes and phases. The curvature of the shocks varies as 
a function of phase due to the elliptical orbit and is also 
different for the two stars. The curvature is stronger 
around the secondary star, with the lower mass-loss rate 
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Fig. 12. — Density in large-box simulation with unequal mass- 
loss rates on elliptical orbits (a) in the orbital plane, (6) along the 
semi-major axis normal to the orbital plane, and (c) normal to the 
previous two slices, at phase 0.5. The same slices, (d), (e), and (/), 
respectively, at phase 0.75, and (g), (h), and (i), respectively, at 
phase 1.0. The color scale is logarithmic, with a^urj—ip—ie from 
1 (white) to 700 (red). The tick marks are spaced aAU. 




Fig. 13. — Temperature in large-box simulation with unequal 
mass- loss rates on elliptical orbits (a) in the orbital plane, (b) along 
the semi- major axis normal to the orbital plane, and (c) normal to 
the previous two slices, at phase 0.5. The same slices, (d), (e), and 
(/), respectively, at phase 0.75, and (g), (h), and (i), respectively, 
at phase 1.0. The color scale is logarithmic, with "^Tjjcv 
from 0.05 (white) to 2.5 (red). The tick marks are spaced oAU. 

wind. The post-shock material between the stars is also 
shifted toward the secondary. 

We will now investigate separately the effect on post- 
shock temperatures, free-free emission, and column den- 
sity of elliptical orbits and unequal winds. In H3.41 we 
present the analysis of a simulation with identical winds 
on an elliptical orbit. In §3.51 we analyze two simula- 
tions, one with differing mass-loss rates and the other 
with differing stellar wind velocities, both on circular or- 
bits. These simulations are run in a smaller box, with 
L = 2.5a AU. 

3.4. Identical Winds on Elliptical Orbits 

We consider here stars with identical winds on el- 
liptical orbits, as described in ^2.2\ Since we do not 
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Fig. 14 



Frequency-integrated free-free power emitted from 
gas with T > 1.155/i!^ keV as a function of phase from SE2.5 
(dotted) and i?2.5 (short dash). The time-average for C2.5 (solid), 
SE2.5 (dash dot), and i?2.5 (long dash) are also shown. Varying 
the stellar separation without full orbital motion is not a good 
approximation to an elliptical orbit. 

model the radiative driving of the winds, we find that the 
highest post-shock temperatures are caused by the wind 
that leaves the masks when the star have their greatest 
velocity towards each other. Pittard (1998), however, 
found that the highest post-shock temperatures occur 
near apastron when the winds have had the longest dis- 
tance to accelerate. Using equation we can compare 
the magnitude of these two effects. Assuming the physi- 
cal parameters for the binary 0-star system described in 
ij2.2l and taking the stellar radius to be i?, = 20i?o, we 
find that the wind velocity at the stagnation point when 
the stars are at apastron should be 6.8% higher than 
when the stars are at periastron. We also find that the 
combined radial and wind velocity is 18% higher when 
the stars are moving toward each other than when they 
are moving away from each other with their highest rel- 
ative velocity. The magnitudes of these effects are com- 
parable. 

We use a temperature cut of T > 1.155/2j/^ keV for the 
elliptical-orbit analysis. The maximum in power radiated 
by free-free emission from gas with T > lAdbftv^ keV 
occurs before periastron, at phase 0.44, when full orbital 
motion is allowed, as shown in Figure [TH While the 
mean power emitted differs by only 3.6% between the cir- 
cular and elliptical cases, the instantaneous value varies 
greatly. The maximum is 79% higher and the minimum 
is 50% lower than the mean. When we instead vary the 
stellar separation without full orbital motion {SE2.5), 
we find that the maximum in power emitted occurs after 
periastron. The time delay from periastron is due to the 
wind released by the star taking time to reach the shock 
front. The mean for SE2.5 is 5.3% lower than for E2.5. 
The maximum for SE2.5 is only 27% higher than the 
mean and the minimum is only 17% lower. The curve 
for SE2.5 is symmetrical in time while that for E2.5 is 
not. We will discuss the reason for this below. 

While the stellar separation is the same at phase (p as at 
phase 1.0 — ip, the stellar radial velocities differ in sign, so 
we do not expect the post-shock temperature distribution 
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Fig. 15.— Volume with T > l.lSSAii'^keV as a function of phase 
for SE2.5 (dotted) and £2.5 (short dash). 

to be symmetrical in time. We instead find that, while 
the velocity of the stars towards each other is increasing, 
the high-temperature tail of the distribution drops off 
more steeply than for a circular orbit and, conversely, 
when the stellar radial velocity is increasing, the high- 
temperature tail drops off less steeply. The sign of the 
change in volume at a given temperature switches first 
for low temperatures and later for higher temperatures. 

While the total volume with T > 1.155/2j^^ keV, shown 
in Figure [T51 in SE2.5 is correlated with the separation 
between the stars, with more volume for larger separa- 
tions, in £'2.5 it is correlated with the stellar radial ve- 
locity, with the most volume when the stars are moving 
towards each other with the highest velocities. At apas- 
tron the shock thickness along the line of centers is ap- 
proximately 14% larger than for the circular-orbit case 
and at periastron it is approximately 14% smaller. The 
phase delay for the gas leaving the mask to reach the 
shock in SE2.5 is typically between 0.02 and 0.04. 

Although the maximum and minimum post-shock tem- 
peratures do not occur near periastron and apastron, re- 
spectively, in our simulation, the maximum and mini- 
mum peak power/volume due to free-free emission do. 
The maximum peak power/volume in E2.5 matches that 
in C2.5 near phases 0.27 and 0.71. At phase 0.27, 
the power/ volume distribution is similar for C2.5 and 
E2.5 except at low values, where there is more vol- 
ume from T > 1.155/ii^^ keV gas for E2.5. At phase 
0.71, the power/volume distribution is again similar ex- 
cept at low values, where there is now less volume from 
T > 1.155fliy'^ keV gas for £2.5. 

In Figure [TBI we compare the column density as a func- 
tion of phase for C2.5 and £2.5. The projection of the 
line of sight onto the orbital plane runs parallel to the 
semi-major axis of the elliptical orbit. When the orbit is 
circular, the two peaks in column density when looking 
toward the center of mass of the system are identical, but 
when the orbit is elliptical, the peak before periastron is 
higher and the peak after is lower. The minimum be- 
fore periastron is lower and the minimum after is higher. 
The difference between the peak and minimum column 
density is 60% higher for the elliptical orbit than for the 



Fig. 16. — Column density on the line of sight described in i|2.3l 
to the center of mass (top), primary star (middle), and secondary 
star (bottom) for C2.5 (dotted) and i?2.5 (short dash). 

circular orbit. 

Since the separation between the stars is no longer 
the same half an orbit apart, the column density to 
the primary and secondary no longer share the same 
phase-shifted time dependence even though the individ- 
ual winds are still identical. The column density to the 
primary is now lower at all phases than it was in the 
circular-orbit case, while the column density to the sec- 
ondary is higher for most, but not all, phases. Where the 
column density to the primary is nearly flat under a cir- 
cular orbit, it now varies with time due to the changing 
stellar separation. The secondary has a narrower peak 
near periastron and isn't quite as flat after apastron. The 
difference between the minimum and peak column den- 
sities is 23% higher for the elliptical orbit than for the 
circular orbit. 

3.5. Unequal Winds on Circular Orbits 

When we consider stars with differing wind parameters 
on circular orbits, we find that the thickness of the post- 
shock region along the line of centers does not change. 
The post-shock region is shifted toward the secondary, 
though, giving a distance ratio from the star centers to 
the center of the post-shock region of 1.2 for CM2.5 and 
1.3 for CW2.5. We expect that modifying the wind ve- 
locity will change the temperature distribution, but that 
modifying the mass- loss rate will not. This follows from 
the equation for energy per unit mass. 



Since T cx P/p, this implies that T (x v^. The wind 
density is set by both the mass-loss rate and wind veloc- 
ity, however, so we expect that modifying either one will 
cause a change in the emissivities. 

We histogram the volume in bins of width AT = 
0.0212/xi^^ keV in post-shock temperature using a tem- 
perature cut of r > 0.9625/ii^^ keV for the unequal 
wind analysis. As expected, the temperature distribu- 
tion changes very little when only the mass-loss rate 
is modified. The total post-shock volume above T — 
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0.9625/11^^ keV in CAI2.5 is lower by 1.3% compared to 
C2.5. The difference in volume in each bin is roughly 
constant over 0.9625 < TkovM"^^^"^ < 1-9425. When 
the wind velocity is modified, however, there is much less 
volume at a given temperature above the cutoff, and the 
peak temperature is slightly lower. Since the secondary 
star has a lower wind velocity, the temperature of the sec- 
ondary's shocked wind is much lower than that of the pri- 
mary. Due to the relation above, we expect the shocked 
wind of the secondary to fall below the temperature cut- 
off. CW2.5 should then have the same temperature dis- 
tribution above T = 0.9625/iJ/^ keV as C2.5, but with the 
volume in all bins reduced by a constant factor. Since 
there is 55% less post-shock volume at all temperatures 
above T = 0.9625/2;^^ keV, there is also 55% less volume 
in each bin in the range 0.9625 < TkovM"^^^"^ < 1-9425. 

We histogram the power per unit volume from gas with 
T > 0.9625/Ii^^ keV and find that the volume peaks at 
low power/volume. CM2.5 has less volume per bin than 
C2.5 except at the very lowest and highest bins. With the 
exception of near and above the highest power/volume 
found in C2.5, CW2.5 also has less volume per bin. 
As with the temperature distribution, the reason for 
this is that the secondary's shocked wind is at too low 
a temperature to be included in the analysis, result- 
ing in less emission. CW2.5 has a much higher maxi- 
mum power/volume than C2.5 and CM2.5. The time- 
averaged power radiated by free-free emission from gas 
with T > 0M25flv^ keV is 0.33% lower for CM2.5 and 
64% lower for CW2.5 than for C2.5. The standard devi- 
ation in this value is 0.58% of the mean for C2.5, 0.30% 
for CM2.5, and 2.8% for CM^2.5. The KH instabihty 
appears to be strongest in the case where the wind ve- 
locities are unequal. 

The column density should be different for both the 
modified wind velocity and mass-loss rate due to changes 
in the density. The column density over a full orbital pe- 
riod for these two cases is shown in Figure [TT] The col- 
umn density to the center of mass is not shown here be- 
cause the post-shock region is no longer centered on that 
location. The mass density of the secondary is lower in 
CM2.5 and higher in CW2.5 compared to C2.5, result- 
ing in a lower and higher column density, respectively, to 
the secondary at all phases. When the column density to 
the primary is at its minimum value, it is the same for the 
equal and unequal wind simulations because the line of 
sight passes only through the wind of the primary, which 
is unchanged. When the line of sight passes through 
both winds, the column density is higher in CW2.5 and 
lower in CM2.5. The column density is at its minimum 
value for the shortest amount of time when the winds 
are equal. Since the shock will wrap more tightly around 
the secondary than the primary, the line of sight will pass 
through only the wind of the primary for a larger range 
of phases when the winds are unequal. 

4. DISCUSSION 

WR 20a is a binary system of nearly-identical Wolf- 
Rayet stars, likely type WN6ha, on circular orbits with 
a period of 3.686 days (Rauw et al. 2005). The primary 
has an estimated mass of 82.7 Mq and the secondary 
81.9 Mq, giving a mass ratio of 1.01. This gives a semi- 
major axis for the system of a = 0.256 AU, or a = 0.256. 
If we assume instead a mass ratio of 1.0, taking Mi — 
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Fig. 17. — Column density on the line of sight described in i|2.3l 
to the primary star (top) and secondary star (bottom) for CM2.5 
(dotted), CW2.5 (short dash), and C2.5 (long dash). 

M2 = 82.3Afo, and use Rauw's assumed wind velocity 
at the stagnation point of 500kms~^ for each star, or 
z/ = 0.5, this would give Vw/Vo « 1.3. We also take 
Rauw et al.'s smooth- wind mass-loss rate of Mi — M2 — 
2.5 X 10~^Mq yr~^, or 77 = 2.5. We also assume fl — 1.4 
(WN). 

Although Vw/Vo for this system is smaller than any 
ratio we modeled, we can still put a lower bound on the 
importance of the Coriolis force by comparing models S 
and C2.5 for this system. The wind mach numbers we 
used in our simulations are smaller than the actual values 
for this system, but since we are already in the strong- 
shock limit, this shouldn't have a qualitative effect on our 
results. For C2.5, we find the power radiated by free-free 
emission to be 5.5 x 10^'^ ergs~^ for T > 0.22 keV. For 
S, we find only 4.8 x lO^'^Z^ ergs^i. This is a difference 
of 13%, even though we used a larger Vw/Vo than is ap- 
propriate. Without taking wind acceleration or radiative 
effects into account, we find that the orbital motion has 
a significant, although not easily measured, effect on the 
power radiated by free-free emission. Unfortunately, by 
eq. ()12p . we find arad — 373 ^ a, meaning this system 
is actually quite radiative. 

Had our simulations incorporated radiative driving of 
the wind, resulting in positive net wind acceleration, the 
increased wind velocity would likely cause the thickness 
of the post-shock region to decrease, with the largest 
effect furthest from the line of centers. The modified 
velocity profile would also give an effective Vw/Vo that 
increases with distance from the line of centers, which 
might also cause the shocks to wrap more tightly around 
the stars. The former effect would cause the wind to 
impact the shock front at a more oblique angle, while 
the latter would have the same effect on the trailing side 
of the shock and the opposite effect on the leading side. 
These two effects, combined with the increased wind ve- 
locity, makes determining the change in the post-shock 
temperature distribution difficult. 

For an isolated star, gravity acts against the radiative 
driving of the wind, leading to a different velocity pro- 
file in the unshocked wind, with lower velocities than 
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those expected without gravity. When these modified 
winds colHde, the thickness of the post-shock region far 
from the line of centers would be increased significantly, 
leading to a larger post-shock volume. Gravity acts not 
only on the unshocked winds, but also on the post-shock 
gas, however. In the post-shock region between a pair 
of equal-mass stars that are held stationary, the force 
due to gravity would pull the post-shock gas toward the 
stagnation point. This would lead to an enhanced den- 
sity near the plane normal to the line of centers, which 
would likely result in a decrease in temperature near the 
contact surface. This discussion is based on simple expec- 
tations; a complete exploration of the effect of gravity on 
the post-shock flow will require 3D simulations in which 
it is included. 

Since a higher ratio of V^/Vo in our simulations cor- 
responds to an increased stellar mass, the effects due to 
gravity will increase with decreasing orbital period. We 
find that the escape velocity for Ku/V"o = 2.5 at a dis- 
tance of 1.25a AU from the stagnation point is roughly 
equal to the unshocked wind velocity at the stagnation 
point. Clearly gravity will have a strong influence on the 
post-shock region in the absence of radiative driving for 
systems with a large Coriolis force, therefore conducting 
simulations that include both the effects of gravity and 
radiative driving is the logical next step. This, however, 
is beyond the scope of this paper. 



Our analysis in ij3.2l shows that, although the orbital 
motion has little effect on the peak temperature and 
density along the line of centers between the stars, it 
can have a significant effect on the temperature and 
emissivity of the gas elsewhere in the post-shock region. 
It also impacts the light curves by changing the time- 
dependence of the attenuation and allows for the K-H 
instability even when the winds are identical. The shape 
of the shocked region changes substantially when orbital 
motion is added, causing the distance between the stars 
and the shock fronts to decrease. The shocks also have 
a rotational velocity once orbital motion is considered, 
causing a difference in the shock strength on the leading 
and trailing sides. While we did not consider radiative 
effects or gravitational forces that would cause the wind 
velocity to vary with distance from the star, adding these 
effects will likely have a large impact on the system. It 
therefore would be fruitful in the future to consider fully 
three-dimensional hydrodynamical models which include 
both radiative driving and gravity. 
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APPENDIX 
WIND GENERATION METHOD 

The wind is generated by imposing the wind solution, described in §2.11 onto a sphere of radius tq, called the masked 
region, before every time step. The values of the density, pressure, and velocity vector within each cell in the masked 
region are computed using the volume-average of the analytic wind solution over the cell, for example for the density 

ip) - (Al) 

where dV = dxdydz is the volume of the cell and r"^ = + + . A similar formula is used for the pressure. The 
radial unit vector, r, is also averaged in this manner, and is used as the direction of the wind velocity in each grid cell 
within the mask in the reference frame of the star. Since the stars are moving, in the reference frame of the grid the 
wind velocity is then 

V = w„(r) -l-v,, (A2) 

where is the stellar velocity. The momentum density vector and the energy density are then calculated from the 
cell-averaged pressure, mass density, and velocity using 

P = (P>v (A3) 

and 

E=^ + \{p)^^\ (A4) 
7 — 1 2 

The singularity in the density and pressure at the center of the star and the issue of dynamic range in the simulation 
are handled simultaneously by imposing a maximum density and pressure, Pmax = 4po a-nd Pmax — (2)^*^/^/0: before 
averaging. These maxima fall at r = O.Srg. 

The integrals of the form of eq. (|Aip are approximated using quadratures by dividing each cell into 10'^ segments 
and then using the midpoint rule. Since space is discretized on the grid, cells which fall only partially within the 
masked region are treated as such. The new grid values for these cells are a linear combination of the previous value 
and the masked value, 

(fn+l = fqmask + (1 " (A5) 

where q = {p, i?,pi,p2,P3) and / is the fraction of the cell's volume covered by the mask. The masks are moved over 
the grid in the appropriate orbit using an N-body integrator. 
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